Superconvergence of Finite Element Discretization of Time Relaxation Models of Advection

The nodal accuracy of …nite element discretizations of advection equations including a time relaxation term is studied. Worst case error estimates have been proven for this combination by energy methods. By considering the Cauchy problem with uniform meshes, precise Fourier analysis of the error is possible. This analysis shows (1)the worst case upper bounds are sharp, (2)time relaxation stabilization does not degrade superconvergence of the usual FEM, and (3)higher order time relaxation is preferable to maintain small numerical errors. Key words. superconvergence, time relaxation, deconvolution 1. Introduction. We consider an approach to eliminating oscillations (forcing them to decay rapidly in time) induced by unresolved scales in conservation laws and convection dominated problems. To reduce the problem to its simplest form (which permits a more exact analysis) consider the advection equation: …nd u = u(x; t) de…ned for x 2 R; t 0 and satisfying ut = ux; 1 < x <1; 0 < t T; (1.1) u(x; 0) = f(x); 1 < x <1: Let over-bar denote a local averaging over radius O( ) (de…ned precisely in Section 1.2). Thus, given an approximate solution u its average is denoted uh and the ‡uctuation is (uh)0 := u uh. Let S (R) denote a …nite element space of smoothest splines de…ned on a uniform mesh (Section 2). The zeroth order example of the approximations we consider is: given a parameter > 0, …nd u : [0; T ] ! S (R) satisfying (ut ux; v) ((uh)0; v) = 0;8v 2 S (R); (1.2) u(x; 0) = I(f)(x); where I is the usual spline interpolation operator. This is the usual Galerkin approximation plus a time relaxation/stabilization term intended to drive small ‡uctuations to zero exponentially fast, see section 1.1. The variational multi-scale framework (see Hughes, Mazzei and Jansen [HMJ00]) gives some insight into this mechanism. Brie‡y, let r(w) := wt wx denote the residual of w in (1.1) and decompose u as u = uh + (uh)0. Setting alternately v = vh and v = (vh)0 in (1.2) gives the equivalent coupled system (ut ux; vh) = (r((uh)0); vh) = 0;8vh 2 Sh; (R); ((u)t (u)x; (vh)0) ((uh)0; v) = (r(uh); (vh)0);8(vh)0 2 (S (R))0: The second equation suggests that the relaxation term will tend to derive ‡uctuations to zero while the …rst suggests that its e¤ects on the means will be limited to the projected residual term on the RHS. There are various other realizations of the same idea. For example, when the map u ! (uh)0 is not positive semi-de…nite, (as can arise, e.g., [LS06] and Pruett [P06]), Department of Mathematics, University of Pittsburgh, Pittsburgh, PA 15260, USA, [email protected], www.math.pitt.edu/~wjl . Partially supported by NSF grant DMS 0508260. 1 the relaxation term should be instead ((uh)0; (vh)0) . The most important variant, analyzed herein and introduced by Stolz, Adams and Kleiser in their computations of turbulent compressible ‡ows [AS02], [SA99], [SAK01a], [SAK01b], [SAK02] (see also [Gue04]), is a higher order time relaxation operator. Brie‡y, (see Section 3.2) given a deconvolution operator DN , i.e., a bounded linear operator on L(R) with the property = DN +O( ) for smooth , (1.3) the higher order1 , generalized ‡uctuation is (u) := u DNu: The higher order time relaxation discretization is then: given > 0 , …nd u : [0; T ]! S (R) satisfying (ut ux; v) ((u) ; v) = 0;8v 2 S (R); (1.4) u(x; 0) = I(f)(x); Note that since = +O( ) (1.2) is the N = 0 case of (1.4). If = 0 , (1.4) reduces to the usual FEM which superconverges at the nodes with rate O(h ). We show that the added stabilization term in (1.4) preserves this property. In Theorem 3.10 we prove u(t) u(t) l2;h Cfh jjf jjH2 +1(R)+ [jjf jjH2N+2(R)+h jjf jjH2N+2 (R)]g: 1.1. The genesis of the time relaxation term. The time relaxation term combines a numerical regularization with a physical model. Because of that it is particularly interesting computationally: it induces a deviation from an exact discretization of (1.1) intended to move the computed solution’s behavior closer to the behavior of the physics (1.1) is often intended to model. In theoretical work on the derivation of conservation laws, regularizations of Chapman-Enskog expansions in Rosenau [R89], Schochet and E. Tadmor [ST92] produced conservation laws with a time relaxation term. This added time relaxation operator is a lower order perturbation and thus (since the equation does not change order or type) questions of well-posedness and boundary conditions are transparent.; in combination with a large eddy simulation model, it has produced positive results for the Navier-Stokes equations at high Reynolds numbers and a mathematical foundation for its inclusion in models for turbulent ‡ow has been derived, [LN05], [LN05b], and [ELN06] . It can also be used quite independently of any turbulence model (and has been so used in compressible ‡ow calculations). As a stand alone regularization, it has been successful for the Euler equations for shock-entropy wave interaction and other tests, [AS02], [SAK01a] , [SAK01b], [SAK02], including aerodynamic noise prediction and control, Guenan¤ [Gue04]. It was observed to ensure su¢ cient numerical entropy dissipation for numerical solution of conservation laws, Adams and Stolz [AS02], p.393. 1.2. Averaging by discrete di¤erential …lters. We study herein averaging by a discrete di¤erential …lter2 (Germano [Ger86] and Manica and Kaya-Merdan 1As N increases to moderate values ! DN becomes quite close to sharp spectral cuto¤, see the …gures in [LN05b]. 2The "best" …lter depends on the neds of application at hand. Scale space analysis suggests the gaussian …lter as the generic case. The above di¤erential …lter arises as the …rst subdiagonal Padé approximation to it in wavenumber space, [GL00]. It is also very convenient for both mathematical analysis and FEM implementation. 2 [MM06]). Germano’s proposal of di¤erential …ltering for ‡uid velocities plays a key role in a number of models of turbulence including the alpha-model, [FHT01], the zeroth order model, [LL03], [LL06a], [LL06b], and deconvolution models, [AS02],[SAK01a] , [SAK01b], [SAK02] and [LMNR06]. Let be the user selected averaging radius (typically = O(h) in computations). Given 2 L(R) its discrete average 2 S (R) is the unique solution of ( x; v h x) + ( ; v ) = ( ; v);8v 2 S (R): (1.5) Associated with (1.5) de…ne the discrete Laplacian operator 4 : L(R) ! S (R) and projection operator h : L(R)! S (R) by ( x; v h x) = ( 4 ; v);8v 2 S (R) , and ( ; v) = ( h ; v);8v 2 S (R): With these de…nitions, the discrete …lter (1.3) can be written = ( 4 + ) ( h ) , or (1.6) ( 4 + ) = ( h ): (1.7) 2. Notation and preliminaries. The fundamental connection between the Galerkin method with splines and an associated di¤erence scheme at the nodes was made by V. Thomée in the early 1970’s in [T72], [T73], see also [TW74], papers of mathematical power and beauty. We shall use the techniques introduced in these papers and shall thus follow the notation in them closely. De…ne, following, e.g., Schoenberg [Sc73], Thomée [T73], the B-spline of order 2. Let denote the characteristic function of [ 1; 1]. De…ne and l by = and l (x) = (h x l). We take S (R) to be the space of splines of at most power growth: S (R) = f X l cl h l (x) : cl = O(jlj) as jlj ! 1 for some qg: The splines in S (R) are C 2 functions which reduce to polynomials of degree 1 on each interval [jh; (j + 1)h] for even and on [(j 1 2 )h; (j + 1 2 )h] for odd. The usual spline interpolation operator is denoted Ih, that is, Ih(v) is that element of S (R) satisfying Ih(v)(lh) = v(lh); l 2 Z: For 2 and integer , 0 2 2, de…ne the trigonometric polynomials (scaled to be independent of h) g ; ( ) = h ( i) 2 1 X l= 1 ( @ @x 0 ; @ @x l )e il , where = [ ]. For even, g ; ( ) is a real, positive trigonometric polynomial, [T73]. Norms associated with doubly in…nite sequences will be useful. For h > 0, the l2;h norm of a function v(x) and sequence c = (cj)j= 1 are de…ned by jjvjj2l2;h = h X j2Z jv(jh)j and jjcjj2l2 = X j2Z jcj j: For c 2 l2 the discrete Fourier transform of c is ec( ) = (Fc)( ) =Pj2Z cje . 3 3. Superconvergence at the nodes. First consider (1.2), i.e., the case N = 0 in (1.4). Expand u(x; t) = P j2Z cj(t) h j (x). Taking the discrete Fourier transform (denoted by an over-tilda) of (1.2) gives hg ;0(h ) d dt ec( ; t) ig ;1(h )ec( ; t) + [hg ;0(h )ec( ; t) ^ (uh; )( ; t)] = 0: (3.1) The di¤erence between various methods of …ltering/averaging and between the time relaxation discretization and the usual Galerkin method lies in the last term on the LHS of the above. The analysis of this last term will be more compact and clear by identifying two trigonometric polynomials that occur frequently. Accordingly, de…ne, suppressing dependence on all parameters except , !( ) := g ;0( ) ( h ) 2g ;2( ) + g ;0( ) ; d( ) := 1 !( ) = ( h ) g ;2( ) ( h ) 2g ;2( ) + g ;0( ) 3.1. Fourier analysis of discrete di¤erential …lters. Consider …rst the discrete di¤erential …lter. If we write u(x; t) = X j2Z cj(t) h j (x) , u (x) = X j2Z cj(t) h j (x) and (uh)0(x) = X